# Ilk, O. (2021) R Yazilimina Giris, Detay Yayincilik, 3. basim

3*(5+7)    

a=3*(5+7)    
a
(a/18)^3  

1:10
seq(1,10)
seq(1,10,3)

rep(seq(1,10,3),4)
rep(seq(1,10,3),each=4)

a=c(4.523,3.487,5.629,2.711)   
a
log(a) 
log10(a)
log(a,10)
1/a
mean(a)  
sd(a)		
sqrt(var(a))	
max(a)
min(a)
which.max(a)
which.min(a)
range(a) 
rev(a) 

a=matrix(c(1.2,3.5,4.7,1.8,-6.2,5.3),3,2)  
a
a=matrix(c(1.2,3.5,4.7,1.8,-6.2,5.3),ncol=2)   
a
a[1,2]
a[3,]    			
a[,2]				
dim(a)			
dim(a)[1]	
nrow(a) 		
dim(a)[2]
ncol(a)			
length(1:10)			

t(a)				
b=t(a)%*%a 		
b
rbind(b,a)				
cbind(b,a) 			
cbind(b,t(a))

d=solve(b)			
d
b%*%d			
library(MASS)
ginv(b)			
eigen(b)			
eigen(b)$values
eigen(b)$vectors
eigen(b)$va
eigen(b)$ve
eigen(b)$v
det(b)  
diag(b)
diag(3) 
diag(5,3,3) 
scale(a,center=T,scale=F) 
mean(a[,1])	
mean(a[,2]) 		
apply(a,2,mean)

obje1=1:3
obje2=4:7
obje3=list(obje1,obje2)
obje3
sapply(obje3,function(x) mean(x))
lapply(obje3,function(x) mean(x))

obje1=1:10
obje2=rep(c("+","-"),each=5)
obje3=data.frame(obje1,obje2)
obje3
tapply(obje3$obje1, obje3$obje2, function(x) c(mean(x), sd(x)))
         
scale(a) 

apply(a,1,sum)		

apply(a, 2, sort)		 

d=matrix(c(1:4,0,6,8,5,3,9,5:9,3:7),ncol=2,byrow=T)
d
d=cbind(1:10,d)
d

d[order(d[,2]),]  

a=c(rep(1,10),100)
a
mean(a)
mean(a,trim=0.1)

ay=factor(1:12, labels="Ay")
ay
gunler=factor(c("pzt","sali","crs","prs","cuma","cmts","pz"))
gunler

cinsiyet=gl(2, 3,6,labels=c("Kadin","Erkek"))
cinsiyet
cinsiyet=gl(2, 3,labels=c("Kadin","Erkek"))
cinsiyet

which(cinsiyet=="Kadin")
class(ay)
ay=data.frame(ay)
class(ay)

x=4:7
y=2:5
z=matrix(c(x,y),ncol=2)
z
class(z)
z[,1]
z$x

w=data.frame(x,y)
w
w[,1]
w$x

for(i in 1:5) print(i)

i=1
repeat{
print(i)
i=i+1
if(i>5) break
}

i=1
while(i<=5){
print(i)
i=i+1} 

for(i in 1:5){
j=1
while(j<3){
ifelse(i==j,print(c(i,j,"Esit")),print(c(i,j,"Esit degil")))
if(i==j) break else j=j+1
}
}

a=read.table("veri1.txt")           
a=matrix(scan("veri1.txt"),ncol=2,byrow=T) 
a

a=read.csv("veri1.csv",header=F) 

a=read.csv(file.choose()) 
a=read.table(file.choose())

a=read.table(file="clipboard",header=F) 

library(foreign)
??foreign::read
?read.spss

write.table(a,"veri2.txt",sep="\t")

write.table(a,"veri2.txt",append=T,sep="\t",row.names=F,col.names=F)

for(i in 1:20){
k=2*i
write.table(k,"cikti1.txt",append=T,sep="\t",row.names=F,col.names=F)}

for(i in 1:20){
k=2*i
write.table(k,"cikti2.txt",sep="\t",row.names=F,col.names=F)}

getwd()  
setwd("c:/")

set.seed(123) 
a=runif(100)		
a=rnorm(100,2,3) 		
a=rpois(100,5)			
a=rbinom(100,10,0.1)		

set.seed(123)
a=rbinom(100,10,0.5) 
length(a) 
a 
a[2] 
a[1:5] 
a[c(1,4,7)] 
b=a[-2] 
length(b) 
b=a[-(1:5)] 
length(b) 
b=a[a>5] 
b
b=a[a==5|a==2] 
b
b=a[a>0&a<4] 
b

index=1:length(a) 
d=cbind(index,a) 
d
b=d[d[,2]>5,] 
b

b=subset(d,a>5) 

sample(1:20,5,replace=T) 

split(1:10,rep(c("Kontrol","Tedavi")),5)  

split(sample(1:10),rep(c("Kontrol","Tedavi")),5) 

x=rnorm(10,-1,3)
x
x[x<0]=0
x

a=rbinom(20,4,0.5)
a
for(i in 1:length(a)) a[i]=a[i+1]
a

i=1:20
a[i]=a[i+1]
a

x=c("abc","abc","dc","ef","kl")
x
xtekil=x[!duplicated(x)]
xtekil
xtekil=unique(x) 
xtekil

install.packages("tidyverse")
library(tidyverse)

set.seed(123)
a=rbinom(5,10,0.1) 
a2=rnorm(5)
index=1:length(a) 
d=data.frame(index,a,a2) 
d

filter(d,a==2) 
arrange(d,a) 

select(d,index,a) 
select(d,-a2) 
select(d,starts_with("a")) 

mutate(d,a3=a**2) 

grup=group_by(d,a) 
summarise(grup,a.ya.gore.a2nin.ortalamasi=mean(a2)) 

e=mutate(d,a3=a**2,a4=c(rep("Ocak.1",each=2),"Ocak.15",rep("Subat.1",each=2))) 
e

ayrik_e <- separate(e,a4, c("Ay", "Gun"))
ayrik_e

ayrik_e <- e %>% separate(a4, c("Ay", "Gun"))
ayrik_e

ayrik_e <- e %>% separate(a4, c("Ay", "Gun"),remove=F)
ayrik_e


x=rbinom(100,10,0.5)
y=rnorm(100,0,1) 
plot(x,y)
plot(x,y,xlab="Binomiyal ölçümler",ylab="Normal ölçümler")
text(3,2.5,"Bu bir metindir")
abline(h=0)

plot(x,y,xlab="Binomiyal ölçümler",ylab="Normal ölçümler",xlim=c(0,9))
abline(v=0)
abline(v=2)

plot(x,y,xlab="Binomiyal ölçümler",ylab="Normal ölçümler")
label=paste("y’nin ortalamasi=",as.character(round(mean(y),2)))
text(5,2.5,label)
text(locator(1), "Bu da baska bir metindir") 
locator(1) 

devAskNewPage(ask = T)
plot(x,y,xlab="Binomiyal ölçümler",ylab="Normal ölçümler")
plot(x,y,xlab="Binomiyal ölçümler",ylab="Normal ölçümler",xlim=c(0,9))

sise=read.table("sise.txt")
names(sise)=c('Sayi','Makina','Teknisyen','Gun')   
boxplot(sise$Sayi~sise$Makina)
boxplot(sise$Sayi~sise$Makina,horizontal=T)


qqnorm(y)
qqline(y)

stem(y)

par(mfrow=c(2,1))
hist(x) 
hist(x,prob=T)


cocuk=read.table("cocuk.txt")
cocuk.table=table(cocuk)
cocuk.table
barplot(cocuk.table,ylab="Frekans",main="Cocuk sayisi") 
barplot(cocuk.table/dim(cocuk)[1],ylab="Olasilik") 
pie(cocuk.table)

icki=matrix(c(41,26,28,9,21,4,9,5),nr=2)
dimnames(icki)=list(Cinsiyet=c("E","K"),Yas=c("15-24","25-44","45-64","65+"))
icki
main=c("Yeni Zelanda'daki Tehlikeli", "Icki Egilimi")
barplot(icki,legend=rownames(icki),ylim=c(0,70),las=1,
main=main,xlab="Yas Grubu",ylab="Frekans")

barplot(icki,legend=rownames(icki),angle=c(60,120),density=10,col="black",
ylim=c(0,70),las=1,main=main,xlab="Yas Grubu",ylab="Frekans")


barplot(icki,legend=rownames(icki),beside=T,las=1,
main=main,xlab="Yas Grubu",ylab="Frekans")

data() 
data(package = .packages(all.available = TRUE)) 

data(faithful)  
attach(faithful) 
hist(eruptions,breaks=15,prob=T) 
lines(density(eruptions)) 

x=seq(10,100,10) 
y1=round(seq(580,450,,10)) 
y2=y1/600 

par(mai=c(1,1,1,1)) 
plot(x,y1,axes=F,xlab="Basinc",ylab="Saglam kalan parca sayisi") 
axis(1,at=seq(10,100,10))
axis(2,at=seq(min(y1),max(y1),10)) 
axis(4,at=y1,labels=round(y2,2)*100) 
mtext("%",side=4,line=2)

x=rnorm(26)
y=2*x+rnorm(26,0,1) 
y[26]=10
veri=data.frame(x,y)
rownames(veri)=LETTERS[1:26]
plot(x,y)
text(x,y, row.names(veri), cex=0.7, pos=2, col="red")

install.packages("plotly")
library(plotly)
plot_ly(x = ~rnorm(100), type = "histogram")

data(iris)
sekil23=plot_ly(iris, x = ~Sepal.Width, y = ~Sepal.Length,type = 'scatter', mode = 'markers', color = ~Species)
sekil23=sekil23 %>% layout(
	xaxis = list(title='Yaprak Genisligi'),
	yaxis = list(title='Yaprak Uzunlugu'))
sekil23


library(ggplot2)

ggplot(iris,aes(x=Species,y=Sepal.Width))+geom_boxplot()+labs(title="Iris verisi için kutu grafigi",
caption = "Veri kaynagı: R demo verisi",
x="Cesit",y="Yaprak Genisligi")

sekil25=ggplot(iris,aes(x = Sepal.Width,y=Sepal.Length, color = Species)) + geom_point()
sekil25=sekil25+labs(x="Yaprak Genisligi",y="Yaprak Uzunlugu")+scale_colour_discrete("Cesit")
sekil25

ggplotly(sekil25)


sqr=function(x) x*x 
sqr(2)

larger=function(x,y){
ifelse(y>x,y,x)	 
} 
larger(1,8)

fact = function(n) if(n<=1) 1 else n*fact(n-1)
fact(4)
factorial
factorial(4)


beta(4,1)
choose(4,1)

harfler=function(not){
if(not<50) "FF"
else if(not<60) "FD"
else if(not<64) "DD"
else if(not<69)  "DC"
else if(not<74)  "CC"
else if(not<79)  "CB"
else if(not<84)  "BB"
else if(not<89)  "BA"
else "AA"
}

harfler(55)
harfler(72)
harfler(95)

newton=function(fun,derf,x0,eps){
iter=0
repeat{
iter=iter+1
x2=x0-fun(x0)/derf(x0)
if(abs(x0-x2)<eps)
break
x0=x2
cat("Iterasyon:",iter,"Sonuc=",x2,fill=T)
}
x2
}

f1=function(x) x^3+sqrt(x)-1
df1=function(x) 3*x^2+(1/2)*x^(-1/2)

x0=seq(0,2,,100)
x0

plot(x0,f1(x0),type="n")
lines(x0,f1(x0))
abline(h=0,lty=2)

newton(f1,df1,0.5,10^-10)

secant=function(fun,x0,x1,eps){
iter=0
repeat{
iter=iter+1
f0=fun(x0)
f1=fun(x1)
x2=(x0*f1-x1*f0)/(f1-f0)
if(abs(x2-x1)<eps)
break
x0=x1
x1=x2
cat("Iterasyon:",iter,"Sonuc=",x2,fill=T)
}
x2
}

f1=function(x) x^3+sqrt(x)-1

secant(f1,0.5,0.55,10^-10)

f2=function(x) x^3-sinh(x)+4*x^2+6*x+9

x0=seq(-10,10,,100)
plot(x0,f2(x0),type="n")
lines(x0,f2(x0))

x0=seq(6,8,,100)
plot(x0,f2(x0),type="n")
lines(x0,f2(x0))
abline(h=0,lty=2)

secant(f2,7,7.2,10^-10)

rlaplace=function(n,beta){
if(beta<=0) print("Beta pozitif bir deger almali") 
k=0
x=NULL
while(k<n){
if(beta<=0) break
k=k+1
u=runif(1)
x[k]=ifelse(u<0.5,beta*log(2*u),-beta*log(2*(1-u)))
}
if(beta>0) x
}

rlaplace(10,2)
rlaplace(10,-2)

rlaplace=function(n,beta){
if(beta<=0) print("Beta pozitif bir deger almali") 
u=runif(n)
x=ifelse(u<0.5,beta*log(2*u),-beta*log(2*(1-u)))
if(beta>0) x
}

rlaplace(10,2)
rlaplace(10,-2)

rucgen=function(n,a,b){
if(b<=a) print("b degeri a degerinden buyuk olmalidir")  
u=runif(n)
x=ifelse(u<0.5,2*a+(b-a)*sqrt(2*u),2*b-(b-a)*sqrt(2*(1-u)))
if(b>a) x
}

rucgen(5,1,5)

set.seed(123) 
a=c(rnorm(10),NA,NA) 
a
mean(a) 
mean(a,na.rm=T) 

a_kayip=is.na(a) 
a_kayip 
sum(!a_kayip)

a_tam=a[!is.na(a)]
a_tam

install.packages('Hmisc')
library(Hmisc)

a_tamamlanan=impute(a) 
a_tamamlanan
median(a_tam)

a_tamamlanan=impute(a,mean) 
a_tamamlanan  

a_tamamlanan=impute(a,0.5)  
a_tamamlanan 

a_tamamlanan=impute(a,"random")  
a_tamamlanan


veri=read.table("satis.txt")
veri
satis=t(veri)
colnames(satis)=c("Yil","Satis","Fiyat","Reklam")
satis

class(satis)

satis=data.frame(satis)
satis
cor(satis) 
shapiro.test(satis[,2]) 
shapiro.test(satis$Satis)

reg1=lm(Satis~Fiyat+Reklam,data=satis) 

summary(reg1) 


names(reg1) 
reg1$fit  

install.packages("car")
library(car)
vif(reg1)


satis.std=satis
for(i in 1:dim(satis)[2]){
satis.std[,i]=(satis[,i]-mean(satis[,i]))/sd(satis[,i])
}

satis.std


satis.std=scale(satis,center=T,scale=T) 
satis.std
class(satis.std)
satis.std=data.frame(satis.std) 



reg2=lm(Satis~Fiyat+Reklam,data=satis.std) 

summary(reg2) 

reg3= lm(Satis~Fiyat-1,data=satis.std) 
summary(reg3) 

anova(reg2,reg3)


# lojistik regresyon
tekrar=read.table("tekrar.txt")
colnames(tekrar)=c("Yas","Zaman","Tekrar")
tekrar
class(tekrar)
cor(tekrar,method="spearman")
reg1=glm(Tekrar~Yas+Zaman,family=binomial,data=tekrar)
summary(reg1)

library(boot)
glm.diag(reg1)
glm.diag.plots(reg1)  
glm.diag.plots(reg1,iden=T) 

cl=ifelse(reg1$fit>0.5,1,0)  
sum(abs(cl-reg1$y)) 
1-sum(abs(cl-reg1$y))/length(reg1$y) 
tablo=table(cl,reg1$y)

install.packages("GMDH2")
library(GMDH2)
confMat(tablo, positive = "1")

library(MASS)
reg1.updated=stepAIC(reg1,trace=T)
summary(reg1.updated)


x=c(5.8,15,21.5,27.5,33.5,39.5,46,51.5)
y=c(0,1,3,8,9,8,10,5)
n=c(98,54,43,48,51,38,28,11)
y/n
sum(y)
sum(n)

grup_lojistik=glm(cbind(y,n-y)~x,family=binomial(logit))
# cbind(sucess,failure)
summary(grup_lojistik)

fitted(grup_lojistik)


reg2=glm(Tekrar~Yas+Zaman,family=binomial(link=probit),data=tekrar)
summary(reg2)

cl=ifelse(reg2$fit>0.5,1,0)  
sum(abs(cl-reg2$y)) 
1-sum(abs(cl-reg2$y))/length(reg2$y) 


sure=c(3.8, 5.3, 3.5, 4.5, 7.2, 5.1)
mean(sure)  
sd(sure) 
length(sure) 

shapiro.test(sure)

t.test(sure,alternative="greater",mu=4)


n=100
x=32

prop.test(x,n,0.5,alternative = "less", correct=F)

makina=read.table("makinalar.txt",header=T)

makina=read.table("clipboard",header=T)
makina

colMeans(makina)
sd(makina$Yeni)


sd(makina$Eski)

shapiro.test(makina$Yeni)
shapiro.test(makina$Eski)

var.test(makina$Yeni,makina$Eski)

# H1: mu_Yeni < mu_Eski
t.test(makina$Yeni,makina$Eski,alternative="less",var.equal=T)

t.test(makina$Yeni,makina$Eski,alternative="two.sided",var.equal=T,conf=0.99)



cinko=read.table("cinko.txt",header=T)
cinko

fark=cinko$alt-cinko$ust


shapiro.test(fark)


t.test(cinko$alt,cinko$ust,alternative="greater",paired=T,conf=0.90)


irs=read.table("irs.txt",header=F)
irs
class(irs)

t.test(irs$V1,mu=160,alt="greater")

qqnorm(irs$V1)
qqline(irs$V1)

shapiro.test(irs$V1)

install.packages("AID")
library(AID)
boxcoxnc(irs$V1)

shapiro.test(irs$V1**0.36)

t.test(irs$V1**0.36,mu=160**0.36,alt="greater")

median(irs$V1)

sum(irs$V1>160)

install.packages("BSDA")
library("BSDA")
SIGN.test(irs$V1,md=160,alt="greater")


bocek_veri=read.table("bocek.csv",header=T,sep=",")
head(bocek_veri)

boxplot(bocek~ilac,data=bocek_veri)

bocek.C=bocek_veri$bocek[bocek_veri$ilac=="C"]
bocek.D=bocek_veri$bocek[bocek_veri$ilac=="D"]

median(bocek.C)
median(bocek.D)

shapiro.test(bocek.C)
shapiro.test(bocek.D)

wilcox.test(bocek.C,bocek.D)


library(MASS)         
head(immer) 

boxplot(immer$Y1, immer$Y2)

wilcox.test(immer$Y1, immer$Y2, paired=TRUE) 

wilcox.test(immer$Y1, immer$Y2, paired=TRUE,alt="greater") 



library(AID)
data(AADT)
?AADT

install.packages("onewaytests")
library(onewaytests)

describe(aadt ~ class, data = AADT)

nor.test(aadt ~ class, data = AADT)

library(MASS)
boxcox(AADT$aadt ~ AADT$class)

AADT$taadt=AADT$aadt**0.1

nor.test( taadt ~ class, data = AADT)

AADT$taadt2=log(AADT$aadt)

nor.test( taadt2 ~ class, data = AADT)

homog.test(taadt2 ~ class, data = AADT)

out=kw.test(aadt ~ class, data = AADT)
paircomp(out)



va1=read.table("VAT1.txt",header=T)
va1

aov1 = aov(Satis~Tasarim,data=va1)  
summary(aov1)
class(va1[,2])  

aov1 = aov(Satis~factor(Tasarim),data=va1)  
summary(aov1)

va1[,2]=as.factor(va1[,2])
aov1 = aov(Satis~Tasarim,data=va1)  
summary(aov1)

install.packages('multcomp')
library(multcomp)
c1=glht(aov1, linfct = mcp(Tasarim = "Tukey"))
summary(c1)

TukeyHSD(aov1, "Tasarim", conf.level=0.9)

library(onewaytests)
aov1 = aov.test(Satis~Tasarim,data=va1)  
# veya aov1 = aov.test(Satis~factor(Tasarim),data=va1)  
paircomp(aov1)



va2=read.table("VAT2.txt",header=T)
va2


interaction.plot(va2[,2],va2[,3],va2[,4], lty=c(1, 12),
                ylab="satis", xlab="uzunluk", lwd=3, trace.label="genislik")
aov2 = aov(satis~uzunluk+genislik,data=va2)  
summary(aov2)
aov2 = aov(satis~uzunluk*genislik,data=va2)  
summary(aov2)
print(model.tables(aov2,"means"),digits=3)

install.packages("twowaytests")
library(twowaytests)
data(alveolar)
descTwoWay(count ~ ovalbumin*treatment, data = alveolar)
nortestTwoWay(count ~ ovalbumin*treatment, data = alveolar)
homogtestTwoWay(count ~ ovalbumin*treatment, data = alveolar, method="Fligner",alpha=0.01)
gpTwoWay(count ~ ovalbumin*treatment, data = alveolar, method = "gPB",alpha=0.01)


va3=read.table("VAT3.txt",header=T)
va3

aov3 = aov(Hatirlama~Etki+Error(Sahis/Etki),va3)
summary(aov3)  

va4=read.table("VAT4.txt",header=T)
va4

aov4=aov(Hatirlama~(Gorev*Etki)+Error(Sahis/(Gorev*Etki)),va4)
summary(aov4)  

va5=read.table("VAT5.txt",header=T)
va5

aov5=aov(Hatirlama~(Gorev*Etki*Cinsiyet*Doz)+Error(Sahis/(Gorev*Etki))+(Cinsiyet*Doz),va5)  
summary(aov5)
a=unlist(summary(aov5)[[4]]) 
length(a)
a[1:5]
a[21:25]

sise=read.table("sise.txt")
names(sise)=c('Sayi','Makina','Teknisyen','Gun')   
sise
class(sise$Makina)
sise$Makina=factor(sise$Makina)
class(sise$Makina)
sise$Teknisyen=factor(sise$Teknisyen)
sise$Gun=factor(sise$Gun)

attach(sise)
aov1=aov(Sayi~ Makina + Makina%in%Teknisyen)
summary(aov1)

aov1=aov(Sayi~ Makina + Makina/Teknisyen)
summary(aov1)
plot(aov1)
detach(sise)

library(MASS)
npk
npk.aov = aov(yield ~ block + N*P*K, npk) 
summary(npk.aov)
summary(npk.aov)[[1]][1,5]
summary(npk.aov)[[1]][,5]  
npk.aov2 = aov(yield ~ block + N+ K, npk) 
summary(npk.aov2)
summary.lm(npk.aov2)  


a=c("-","+")
tasarim1=expand.grid(A=a,B=a,C=a)
tasarim1 
class(tasarim1[,1])

a=c(0,1)
tasarim1=expand.grid(A=a,B=a,C=a)
tasarim1
class(tasarim1[,1])

a=c("0","1")
tasarim1=expand.grid(A=a,B=a,C=a)
tasarim1
class(tasarim1[,1])  


ci=function(x,n,a){

v1=2*x
v2=2*(n-x+1)
v3=2*(x+1)
v4=2*(n-x)

a2=1-a/2
invb1=1
if(v1>0) invb1=qf(a2,v2,v1)

invb2=1
if(v4>0) invb2=qf(a2,v3,v4)
plkesin=v1/(v1+v2*invb1)
pukesin=v3*invb2/(v4+v3*invb2)

cat("kesin CI=", plkesin,pukesin,sep="\n")

z=-qnorm(1-a2)
p=x/n;
pl=p-z*sqrt(p*(1-p)/n)
pu=p+z*sqrt(p*(1-p)/n)

cat("Yaklasik CI=",pl,pu,sep="\n")
}

ci(0,10,0.05)


ci(0,100,0.05)



for(i in 1:10){
x=rgamma(100,1)
y=rgamma(100,1)
result[i,]=coef(lm(y~x))
}




result=matrix(0,10,ncol=2)
result

for(i in 1:10){
x=rgamma(100,1)
y=rgamma(100,1)
result[i,]=coef(lm(y~x))
}

result




reg=function(m){
count=0
while(count<m){
e=rnorm(50,0,1)
x=runif(50,-2,2)
y=10*x+e
lmbeta0=lm(y~x)$coef[1]
lmbeta1=lm(y~x)$coef[2]
lmc=matrix(c(lmbeta0,lmbeta1),ncol=2)
write.table(lmc,"reg_sim_coeffs.txt",col.names=F,row.names=F,append=T)
count=count+1
}
}

reg(20)
res=read.table("reg_sim_coeffs.txt")
res
mean(res[,1])
mean(res[,2])


(mean(res[,2])-10)**2 


var(res[,2])+(mean(res[,2])-10)**2 



n=10
nsim=10000
theta.hat = rep(0,nsim)

for (i in 1:nsim) {
    x = rcauchy(n)
   theta.hat[i] = median(x)
}

hist(theta.hat, freq = FALSE, breaks = 100)



t_stat_sim=function(N,n,alpha){
tval=qt(1-alpha,n-1)
z=matrix(rnorm(N*n),N,n)
tstat=function(x) mean(x)/sqrt(var(x)/length(x))
ahat=sum(apply(z,1,tstat)>tval)/N
shat=sqrt(ahat*(1-ahat)/(N-1))
write(c(alpha,ahat,shat),file="mc_out.txt",ncol=3,append=T)
}

t_stat_sim(10000,10,0.01)
t_stat_sim(10000,10,0.025)
t_stat_sim(10000,10,0.05)
t_stat_sim(10000,10,0.075)
t_stat_sim(10000,10,0.1)

sonuclar=read.table("mc_out.txt")
alpha=sonuclar[,1]
ahat=sonuclar[,2]
shat=sonuclar[,3]
plot(alpha,ahat,xlab="Alpha",ylab="Error Rate",ylim=c(min(ahat-shat),max(ahat+shat)))
segments(alpha,ahat-shat,alpha,ahat+shat)


bnormal=function(N,rho){
rho2=sqrt(1-rho^2)
k=1	
theta1=0
theta2=0
while(k<=N){
theta1=rnorm(1,rho*theta2,rho2)
theta2=rnorm(1,rho*theta1,rho2)
write.table(matrix(c(k,theta1,theta2),ncol=3),"bivariatenormal.txt",sep="\t",append=T,row.names=F,col.names=F)
k=k+1
}
}

bnormal(1000,0.5)   

gen.theta=read.table("bivariatenormal.txt")

par(mfrow=c(2,1))
plot(gen.theta[,1],gen.theta[,2],xlab="Iteration",ylab="Theta1",type='n')
lines(gen.theta[,1],gen.theta[,2])
plot(gen.theta[,1],gen.theta[,3],xlab="Iteration",ylab="Theta2",type='n')
lines(gen.theta[,1],gen.theta[,3])

var(gen.theta[,2:3])   

apply(gen.theta[,2:3],2,mean)   

par(mfrow=c(2,1))
hist(gen.theta[,2],xlab="Theta1",main="Histogram of Theta1")
hist(gen.theta[,3],xlab="Theta2",main="Histogram of Theta2")